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We use multidimensional diffusion processes to approximate the 
dynamics of a queue served by many parallel servers. The queue is 
served in the first-in-first-out (FIFO) order and the customers wait- 
ing in queue may abandon the system without service. Two diffusion 
models are proposed in this paper. They differ in how the patience 
time distribution is built into them. The first diffusion model uses 
the patience time density at zero and the second one uses the en- 
tire patience time distribution. To analyze these diffusion models, we 
develop a numerical algorithm for computing the stationary distri- 
bution of such a diffusion process. A crucial part of the algorithm is 
to choose an appropriate reference density. Using a conjecture on the 
tail behavior of a limit queue length process, we propose a system- 
atic approach to constructing a reference density. With the proposed 
reference density, the algorithm is shown to converge quickly in nu- 
merical experiments. These experiments also show that the diffusion 
models are good approximations for many-server queues, sometimes 
for queues with as few as twenty servers. 

1. Introduction. The focus of this paper is the numerical analysis of 
multidimensional diffusion processes that approximate the dynamics of a 
queue with many parallel servers. A many-server queue serves as a building 
block modeling operations of a large-scale service system. Such a service 
system may be a call center with hundreds of agents, a hospital depart- 
ment with tens or hundreds of inpatient beds, or a computer cluster with 
many processors. When the customers of a service system are human be- 
ings, some of them may abandon the system before their service begins. 
The phenomenon of customer abandonment is ubiquitous because no one 
would wait for service indefinitely. As argued in Garnett, Mandelbaum and 
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Reiman (2002), one must model customer abandonment explicitly in order 
for an operational model to be relevant for decision making. We model cus- 
tomer abandonment by assigning each customer a patience time. When a 
customer's waiting time for service exceeds his patience time, he abandons 
the queue without service. 

The exact analysis of such a many-server queue has been largely limited 
to an M/M/n+M model (also called an Erlang-A model) that has a Poisson 
arrival process and exponential service and patience time distributions. See, 
e.g., Garnett, Mandelbaum and Reiman (2002). However, as pointed out by 
Brown et al. (2005), the service time distribution in a call center appears to 
follow a log-normal distribution. Such distributions have also been observed 
by Shi et al. (2010) for lengths of stay in a hospital. Moreover, the patience 
time distribution in a call center has been observed to be far from exponential 
by Zeltyn and Mandelbaum (2005). With a general service or patience time 
distribution, there is no finite-dimensional Markovian representation of the 
queue. Except computer simulations, there is no method to exactly analyze 
such a queue either analytically or numerically. To deal with the challenge, 
the following strategies are adopted in this paper for analyzing a many-server 
queue. 

First, the service time distribution is restricted to be phase-type. Since 
phase-type distributions can be used to approximate any positive-valued 
distribution, such a queueing model is still relevant to practical systems. 
We focus on a Gl/Ph/n + GI queue with n identical servers. The first GI 
indicates that the customer interarrival times are independent and identi- 
cally distributed (iid) following a general distribution, the Ph indicates that 
the service times are iid following a phase-type distribution, and the +GI 
indicates that the patience times are iid following a general distribution. Sec- 
ond, we are particularly interested in a queue operating in the Quality- and 
Efficiency-Driven (QED) regime: The queue has a large number of servers 
and the arrival rate is high; the arrival rate and the service capacity are 
approximately balanced so that the mean waiting time is relatively short 
compared with the mean service time. As argued in Garnett, Mandelbaum 
and Reiman (2002), such a system has high server utilization as well as 
short customer waiting times and a small fraction of abandonment. There- 
fore, both quality and efficiency can be achieved in this regime. Third, rather 
than analyzing the many-server queue itself, we propose and analyze diffu- 
sion models that approximate the queue. Two diffusion models are proposed 
in this paper. In each model, a multidimensional diffusion process is used to 
represent the scaled customer numbers among service phases. The difference 
between the two diffusion models lies in how the patience time distribution 
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is built into them. The first diffusion model uses the patience time density at 
zero and the second one uses the entire patience time distribution. In partic- 
ular, the diffusion process in the first model is a multidimensional piecewise 
Ornstein-Uhlenbeck (OU) process. We propose an algorithm in this paper 
to numerically solve the stationary distribution of a diffusion process. The 
computed stationary distribution is used to estimate the performance mea- 
sures of a many-server queue. Numerical examples in Section 6 demonstrate 
that the diffusion models are very accurate in predicting the performance of 
a many-server queue, even if the queue has as few as twenty servers. 

Except for the one-dimensional case, the stationary distribution of a piece- 
wise OU process has no explicit formula. The algorithm proposed in this 
paper is a variant of the one in Dai and Harrison (1992), which computes 
the stationary distribution of a semimartingale reflecting Brownian motion 
(SRBM). As in Dai and Harrison (1992), the starting point of our algorithm 
is the basic adjoint relationship that characterizes the stationary distribution 
of a diffusion process. With an appropriate reference density, the algorithm 
can produce a stationary density that satisfies this relationship. 

We set up a Hilbert space using the reference density. In this space, the 
stationary density is orthogonal to an infinite-dimensional subspace H. A 
finite-dimensional subspace is used to approximate H and a function or- 
thogonal to Hk can be numerically computed by solving a system of finitely 
many linear equations. This function is used to approximate the stationary 
density. There are two sources of error in computing the approximate sta- 
tionary density by our algorithm: approximation error and round-off error. 
The approximation error arises because is an approximation of H. As 
Hk increases to H, the approximation error decreases to zero. The round-off 
error occurs because the solution to the system of linear equations has error 
due to the finite precision of a computer. As increases to H, the dimension 
of the linear system gets higher and the coefficient matrix becomes closer 
to singular. As a consequence, the round-off error increases. The condition 
number of the matrix is used as a proxy for the round-off error. Balancing 
the approximation error and the round-off error is an important issue in our 
algorithm. 

A properly chosen reference density is essential for the convergence of the 
algorithm. By convergence, we mean that the approximation error converges 
to zero as increases to H. More importantly, a "good" reference density 
can make Hk converge to H quickly so that the resulting approximation error 
and round-off error are small simultaneously even though the dimension of 
Hk is moderate. To ensure the convergence of the algorithm, the reference 
density should have a comparable or slower decay rate than the stationary 
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density. Since the stationary density is unknown, we make a conjecture on 
the tail behavior of the limit queue length process of many-server queues with 
customer abandonment. We conjecture that the limit queue length process 
has a Gaussian tail and the tail depends on the service time distribution 
only through its first two moments. This tail is used to construct a product- 
form reference density. With this reference density, the algorithm appears 
to converge quickly, producing stable and accurate results. For comparison 
purposes, we also test the algorithm with a certain "naively" chosen reference 
density in Section 7.1. The algorithm fails to converge with the "naive" 
reference density. The major contributions of this paper are the proposed 
diffusion models and the proposed reference densities that are critical to the 
numerical algorithm for computing the stationary distribution of a diffusion 
model. 

Our diffusion models are obtained by replacing certain scaled renewal 
processes by Brownian motions. The replacement procedure is rooted in the 
many-server heavy traffic limit theorems that are proved in an asymptotic 
regime. The two diffusion model proposed in this paper are motivated by the 
diffusion limits proved in Dai, He and Tezcan (2010) and Reed and Tezcan 
(2009). See Section 4.3 for more details. The theory of diffusion approxi- 
mation for many-server queues can be traced back to the seminal paper by 
Halfin and Whitt (1981), where a diffusion limit was established for GI/M/n 
queues. Garnett, Mandelbaum and Reiman (2002) proved a diffusion limit 
for M/M/n + M queues that allows for customer abandonment, and Whitt 
(2005) generalized the result to G/M/n + M queues. Puhalskii and Reiman 
(2000) established a diffusion limit for Gl/Ph/n queues. Their result was 
extended to G/Ph/n + GI queues with customer abandonment in Dai, He 
and Tezcan (2010). Recently, Reed and Tezcan (2009) proved a diffusion 
limit for GI/M/n + GI queues. In their framework, a refined limit process 
is obtained by scaling the patience time hazard rate function. 

Harrison and Nguyen (1990) derived Brownian models for multiclass open 
queueing networks. Their diffusion models are SRBMs and are rooted in the 
conventional heavy traffic limit theorems that are pioneered in Iglehart and 
Whitt (1970) for serial networks and Reiman (1984) for single-class networks. 
See Williams (1996) for a survey of limit theorems in literature. For a two- 
dimensional SRBM living in a rectangle, Dai and Harrison (1991) proposed 
an algorithm computing its stationary distribution. Dai and Harrison (1992) 
extended the algorithm for an SRBM living in an orthant. To deal with the 
unbounded state space, the notion of a reference density was first introduced 
there. Their finite-dimensional space is constructed via (global) multi- 
nominals of order at most k. With this choice of H^, the algorithm appears 
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numerically unstable occasionally. In such a case, the round-off error may 
dominate the approximation error while the approximation error is still sig- 
nificant. Shen et al. (2002) extended Dai and Harrison (1991) to a hypercube 
state space of an arbitrary dimension. They used a finite element method 
to construct to avoid numerical instability. Their algorithm sometimes 
converges slowly because they did not explore a reference density. A linear 
programming algorithm for computing the stationary distribution of a diffu- 
sion process was proposed in Saure, Glynn and Zeevi (2009). Both SRBMs 
in an orthant and a diffusion approximation of many-server queues with 
two priority classes were investigated in their paper. Like the role of the 
reference density, it appears that the rescaling of variables is essential to the 
convergence of their algorithm. 

The remainder of the paper is organized as follows. General diffusion pro- 
cesses are introduced in Section 2, where the basic adjoint relationship for a 
diffusion process is also presented. In Section 3, we begin with recapitulat- 
ing the generic algorithm of Dai and Harrison (1992), and then propose a 
finite element implementation that follows Shen et al. (2002). Two diffusion 
models for Gl/Ph/n + GI queues are presented in Section 4. In Section 5, 
we discuss how to choose an appropriate reference density exploiting the 
tail behavior of a diffusion process. In Section 6, it is demonstrated via nu- 
merical examples that the diffusion models serve as good approximations of 
many-server queues. Section 7 is dedicated to some implementation issues 
arising from the proposed algorithm. The paper is concluded in Section 8. 
We leave the proofs of Propositions 2 and 3 to the appendix. 

Notation. The symbols N, M, and M + are used to denote the sets of 
positive integers, real numbers, and nonnegative real numbers, respectively. 
For d, m E N, R d denotes the d-dimensional Euclidean space and ]^ dxm 
denotes the space oidy.ni real matrices. We use Cf (M d ) to denote the set of 
real-valued functions on M d that are twice continuously differentiable with 
bounded first and second derivatives. For z, w E M, we set z + = max{z,0}, 
z~ = max{— z,0}, and z A w = min{z,w}. All vectors are envisioned as 
column vectors. For a d-dimensional vector x E we use xj for its jth 
entry and diag(x) for the d x d diagonal matrix with jth diagonal entry Xj. 
For a matrix M, M' denotes its transpose, My denotes its (i, j)th entry, 
and \M\ = (%2ij M^) 1 ^ 2 - We reserve / for the d x d identity matrix, e for 
the d-dimensional vector of ones, and e 3 for the d-dimensional vector with 
its jth entry one and all other entries zero. Given two functions (p and (p 
from N to R, we write fi(n) = 0{ip{n)) as n — > oo if there exists a constant 
k > and some no E N such that |<^(n)| < K|(/j(n)| for all n > tiq. 
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2. Diffusion processes. Let d be a positive integer. This paper focuses 
on a d-dimensional diffusion process X = {X(t) : t > 0}. Let (il, J 7 , F, P) be 
a filtered probability space with filtration F = {Tt ■ t > 0}. We assume that 
X satisfies the following stochastic differential equation 

(2.1) X{t) = X(0)+ I b(X(s))ds+ [ a{X{s))dB(s), 

Jo Jo 

where the drift coefficient b is a function from W 1 to the diffusion co- 
efficient a is a function from R d to R dxm , and B = {B(t) : t > 0} is an 
m-dimensional standard Brownian motion with respect to F. We assume 
that both b and a are Lipschitz continuous, i.e., there exists a constant 
c\ > such that 

(2.2) \b(x) - b(y)\ + \a(x) - a{y)\ < d\x - y\ for all x, y G M d . 

Under condition (2.2), the stochastic differential equation (2.1) has a unique 
strong solution, i.e., there exists a unique process X on (f2, J 7 , F, F) such that 
(a) X is adapted to F, (b) for each sample path u G Q, X(t,u)) is continuous 
in t, and (c) for each t > 0, the stochastic differential equation (2.1) holds 
with probability one. See 0ksendal (2003) for more details. We also assume 
that a is uniformly elliptic, i.e., there exists a constant C2 > such that 

(2.3) y'^{x)y > c 2 y'y for all x, y G M d , 
where 

(2.4) S(ar) = a{x)<r'(x). 

We are interested in the diffusion processes that model the dynamics of a 
queue with many parallel servers. Parallel-server queues will be introduced in 
Section 4. In that section, two diffusion processes will be identified to model 
such a queue and the coefficients b and a will be mapped out explicitly 
in terms of primitive data of the queue. The diffusion models presented in 
Section 4 are rooted in many-server heavy traffic limit theorems proved in 
Dai, He and Tezcan (2010) and Reed and Tezcan (2009). 

A probability distribution tt on M rf is said to be a stationary distribu- 
tion of X if X(t) follows distribution tt for each t > whenever X(0) has 
distribution tt. Condition (2.3) is required to ensure the uniqueness of the 
stationary distribution. See Dieker and Gao (2011) for more details. In this 
paper, we assume that X has a unique stationary distribution tt and tt has a 
density g with respect to the Lebesgue measure on For a general diffusion 
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process, there is no explicit solution for tt. This paper develops a numerical 
algorithm computing tt. As in Dai and Harrison (1992), the starting point 
of the algorithm is the basic adjoint relationship 

(2.5) f Qf{x) vr(dx) = for all / € C 6 2 (R d ), 

JR d 

where Q is the generator of X defined by 
(2-6) 

g f ( x) =EM^+^E s ^f^£ for each / G c ^ 
j=i 3 j=i i=\ j 1 

and £ is the covariance matrix given by (2.4). The following theorem is a 
consequence of Proposition 9.2 in Ethier and Kurtz (1986). 

Theorem 1. Let tt be a probability distribution onM. d that satisfies (2.5). 
Then, tt is a stationary distribution of X. 

In this paper, we conjecture that a stronger version of Theorem 1 is true. 

Conjecture 2. Let tt be a signed measure on M d that satisfies (2.5) 
and ir(R d ) = 1. Then, ir is a nonnegative measure and consequently it is a 
stationary distribution of X . 

Our algorithm is to construct a function g on R rf such that 
(2.7) 

/ g(x) dx = 1 and f Gf(x)g(x) dx = for all / G C%(R d ). 

Assuming that Conjecture 2 is true, g must be the unique stationary density 
of X. As a special case, the nonnegativity of a signed measure tt that satisfies 
(2.5) for a piecewise OU process was proposed as an open problem by Dai and 
Dieker (2010). Piecewise OU processes will be introduced in Section 4.3.1. 

3. A finite element algorithm for stationary distributions. In 

this section, we propose a numerical algorithm computing the stationary 
density g. The basic algorithm follows the one developed in Dai and Harrison 
(1992). The finite element implementation closely follows Shen et al. (2002). 

3.1. A reference density. To compute the stationary density g, we adopt 
a notion called a reference density that was first introduced by Dai and 
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Harrison (1992). A reference density for g is a function r denned from R d to 
K-L such that 



(3.1) / r(x) dx < oo and / q (x)r(x) dx < oo, 

where 

q( x ) = f or eac h X £R d 

r(x) 

is called the ratio function. Such a function r exists because r = g is a 
reference density. The reference density controls the convergence of our al- 
gorithm. We will discuss how to choose a reference density for the diffusion 
models of a many-server queue in Section 5. 

For the rest of Section 3, we assume that a reference density r satisfying 

(3.1) has been determined and remains fixed. In addition, we assume that 

(3.2) / b 2 Ax)r(x) dx < co and / T,%(x)r(x) dx < oo 

for j,£ = l,...,d. Since both b and a are Lipschitz continuous, condition 

(3.2) is satisfied if 

(3.3) / |x| 4 r(x) dx < oo. 

jR d 

Let L 2 (M. d ,r) be the space of all square-integrable functions on with 
respect to the measure that has density r, i.e., 

L 2 {R d , r) = {/ G £(M d ) : f f 2 (x)r(x) dx < oo) 

where B(R d ) is the set of Borel-measurable functions on M. d . Condition (3.1) 
implies that q G L 2 (K d ,r). We define an inner product on L 2 (R d ,r) by 

(/,/> = / f(x)f(x)r(x)dx for /,/ G L 2 (R d ,r). 

The induced norm is given by 

(3.4) ll/H = (Z,/) 1 / 2 for each / G L 2 (R d ,r). 

One can check that L 2 (R d , r) is a Hilbert space and assumption (3.2) ensures 
that Qf G L 2 (R d , r) for all / G C 2 (R d ). In L 2 (M d ,r), the basic adjoint 
relationship in (2.7) is equivalent to 



(3.5) (Of, q}=0 for all / € C\ 
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With a fixed reference density r, we need only compute the ratio function 
q by (3.5). Once q is obtained, we can compute the stationary density via 
g(x) = q(x)r(x) for x G 
Let 

(3.6) H = the closure of {Gf : f G C 6 2 (M d )} 

where the closure is taken in the norm in (3.4). As a subspace of L 2 (R d ,r), 
H is orthogonal to q. Let c be a constant function and c(x) = 1 for all 
x £R d . Clearly, c G L 2 (R d ,r) but c £ H because 

(3.7) (c,q) = [ g(x)dx = l. 

JR d 

Let 

(3.8) c = argmin ||c — /|| 

be the projection of c onto H. Then, c— c must be orthogonal to H. Assuming 
that Conjecture 2 holds and X has a unique stationary density g, one must 
have q = n q (c — c) for some constant k 9 G M. By (3.7), the normalizing 
constant n q satisfies 

K~ 1 = (c, c — c) = (c — c, c — c) + (c, c — c) = ||c — c|| 2 . 

Hence, the ratio function is given by 

(3.9) * = 

\\c — c\\ 

3.2. An approximate stationary density. To compute q by (3.9), we need 
first compute c, the projection of c onto H. The space H is linear and 
infinite-dimensional (i.e., a basis of H contains infinitely many functions). 
In general, solving (3.8) in an infinite-dimensional space is impossible. In 
the algorithm, we use a finite-dimensional subspace to approximate H. 

Suppose that there exists a sequence of finite-dimensional subspaces {Hk : 
k G N} of H such that H k -> H in L 2 (M d , r) as fc ->■ oo. Here, H k ^ H in 
L 2 (K d ,r) means that for each f £ H, there exists a sequence of functions 
{(fk : A; G N} with 99^. G .fffc such that ||</?fe — /|| — >• as k — > 00. Let 

c fc = argmin ||c - /|| 
feH k 

be the projection of c onto i^^. By Proposition 7 of Dai and Harrison (1992), 
we have the following approximation result. 
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Proposition 1. Assume that Conjecture 2 is true. Then, 

\\lk — 1 1 — ^ as k — > oo, 

where q k = (c — c&)/ ||c — Cfc|| 2 . Furthermore, when the reference density r is 
bounded, 

(Sk{ x ) ~ g( x )) 2 dx — )• as k — )• oo, 



/ 



where gk{x) = qk(x)r(x) for each x E M. d . 
As in Dai and Harrison (1992), we choose 

(3.10) H k = {Gf:fe C k } 

for some finite-dimensional space C k . In Section 3.3, we will discuss how to 
construct C k using a finite element method. For notational convenience, we 
omit the subscript k when k is fixed. The finite-dimensional functional space 
is thus denoted by C. Let mc be the dimension of C and {/, : i = 1, . . . , mc} 
be a basis of C. We assume that the family {Gfi : i = 1, . . . , mc} is linearly 
independent in L 2 (M. d ,r). Then, 

mc 

(3.11) Cfc = UiQfi for some Ui G K and i = 1, . . . , mc- 

i=i 

Using the fact (Gfi,c — c~k) = for i = 1, . . . , mc, we obtain a system of 
linear equations 

(3.12) Au = v 
where 

(3.13) An = {Gfi, Gfe), u = (ui, . . . ,Um c )', Vi = (Gfi,c). 

By the linear independence assumption, the mc x mc matrix A is positive 
definite. Thus, u = A~ 1 v is the unique solution to (3.12). Once the vector u is 
obtained, we can compute the projection Ck by (3.11). Finally, the stationary 
density g can be approximated via 

g{x) sa g k {x) = r(x) C ^ ~_ Cfc ^ for each x G R d . 
\\ C ~ Cfc 1 1 
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3.3. A finite element method. In Dai and Harrison (1992), the authors 
employed multinominals of orders up to k to construct the space This 
choice appears to be numerically unstable. The approximation error is signif- 
icant when k is small, say, k < 5. As k increases, the round-off error in solv- 
ing (3.12) increases and ultimately dominates the approximation error. Al- 
though their implementation produces accurate estimates for the stationary 
means of SRBMs, it sometimes produces poor estimates for the stationary 
distributions. In this section, we construct a sequence of spaces {Ck '■ k G N} 
using the finite element method as in Shen et al. (2002). Because the state 
space in Shen et al. (2002) is bounded, neither a reference density nor state 
space truncation is used there. 

The state space of X is unbounded in our setting. It is necessary to 
truncate the state space to apply the finite element method. Let {K^ : k G 
N} be a sequence of compact sets in For each / £ Cj, we assume that 
f(x) = for x G M d \ Kj,. The subscript k is omitted again when it is 
fixed and we use K to denote the compact support of the space C. In our 
implementation, we restrict K to be a cf-dimensional hypercube 

(3.14) K= [-Ci,6] x ••• x [-&,&], 

where both Q and £ 3 - are positive constants for j = 1, . . . , d. 

We partition K into a finite number of subdomains. Such a partition is 
called a mesh and each subdomain is called a finite element. Since K is 
a hypercube, it is natural to use a lattice mesh, where each finite element 
is again a hypercube. In this case, each corner point of a finite element is 
called a node. In dimension j = 1, . . . , d, we divide the interval [— Q, £,•] into 
rij subintervals by partition points 

-Cj = Vi<Vj<"-<Vj 

Then, K is divided into IIj=i n j finite elements. For future reference, we 
label the nodes in the way that node (ii, . . . , id) corresponds to spatial co- 
ordinate (y^ 1 , . . . , y d d ) , and define 

h l j = yj +1 - y j for i = 0, . . . , rij — 1 and j = 1, . . . , d. 

If A denotes such a mesh, we define 

| A| = maxj/i' : I = 0, . . . , rij — 1; j = 1, . . . , d} 

and 

\h% \ 

(3.15) 7/a = max j^- : 4,4 = 0, . . . , rij - 1; ji,j 2 = 1,. . . ,d; j\ ^ J2 j- 
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The finite-dimensional space C is generated using the above mesh. We 
use the cubic Hermite basis functions to construct a basis of C, as in Shen 
et al. (2002). The one-dimensional Hermite basis functions for —1 < z < 1 
are given by 

(3.16) <j>{z) = {\z\ - l) 2 {2\z\ + 1) and t/j(z) = z(\z\ - l) 2 . 
In dimension j = 1, . . . , d and for I = 1, . . . , rij — 1, let 



and 



^Az) = 



ui-l 



if Uj < z < y-j 
otherwise 



ifyj _1 <*<l£ 







if y] < z < y) 
otherwise. 



e+i 



Let x = (xi, . . . , Xd)' be a vector in K. At node (ii, . . . , id), the basis func- 
tions of C are the tensor-product Hermite basis functions 



(3-17) f iu .. 
where Xj is either or 1 and 



x 



3=1 



<t> l Hz) ifxi = o, 



^/ 0) if Xj 



1. 



Therefore, each node has 2 d tensor-product basis functions and the space C 
has a total of 



(3.18) 

basis functions. 



m c = 2 d Y[( nj -l) 

3=1 
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The space C is not a subspace of C 2 (Mr) . For the one-dimensional Hermite 
basis functions in (3.16), the second order derivative of <ft(z) is not defined 
at z = — 1 and 1, and the second order derivative of ip( z ) is not defined 
at z = —1, 0, and 1. As a consequence, there exists f £ C for which Qf 
is not defined on the boundaries of certain finite elements. Because such 
boundaries have Lebesgue measure zero in M. d , for each / G C, we can find 
a sequence of functions {<pi : i € N} in Cf(M d ) such that \\Gipi - Gf\\ -»• 
as i — > oo. Hence, C H still holds for each k. 

For the linear system (3.12) to have a unique solution, the family of func- 
tions 

{Gfh,...,i d ,xu-,xd '■ *j = !> - 1; Xj = 0, 1; j = 1,. . .d} 

must be linearly independent in L 2 (R d ,r). The following proposition pro- 
vides sufficient conditions for the linear independence. Its proof can be found 
in the appendix. 

Proposition 2. Let Q be the generator of X in (2.6) such that condi- 
tions (2.2) and (2.3) holds and all entries of E are continuously differen- 
tiable. Assume that r(x) > for all x € M. d . Then, the family of functions 

{Gfii,...,i d ,xi,...,Xd -ij = rij - l;xj =0,1; j = l,...d} 

is linearly independent in L 2 (lR. d ,r), where /ii,...,i d ,xi ...,xd * s ^ e basis func- 
tion of C given by (3.17). Consequently, the solution to the linear system 
(3.12) is unique. 

Now let us consider a sequence of functional spaces {G\. : k E N}. Let 
Afc be the mesh for constructing GV We assume that the mesh A^+i is a 
refinement of A/%, i.e., a node or an interelement boundary in A^ is also a 
node or an interelement boundary in Afc+i. We further assume that such 
refinements are regular, i.e., for each r]A k defined in (3.15), the set {r]A k '■ 
k £ N} is bounded. The next proposition, along with Proposition 1, justifies 
the proposed algorithm for computing the stationary distribution. We leave 
the proof of Proposition 3 to the appendix, too. 

Proposition 3. Let {A^ : k £ N} be a sequence of lattice meshes such 
that each A/% + i is a refinement of A& and the refinements are regular. Let 
Kk be the d- dimensional finite hypercube that is the domain of A}~, andCk be 
the functional space generated by A^ using the tensor-product Hermite basis 
functions in (3.17). Let H be the infinite- dimensional space in (3.6) and 
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Hk be the finite- dimensional space in (3.10), where the generator Q satisfies 
(2.2) and (3.2). Assume that 

|Afc| ->■ and K k t ask^-oo. 

Then, 

Hk — > H as k —> oo. 

4. Diffusion models for many-server queues. In this section, we 
introduce Gl/Ph/n + GI queues and present two diffusion models for such 
a queue with many servers. The two models differ in how the patience time 
distribution is built into them. The patience time density at zero is used in 
the first model, whereas the entire patience time distribution is used in the 
second model. 

4.1. Gl/Ph/n + GI queues in the QED regime. We focus on a queue 
with many servers working in the QED regime. The QED regime will be 
discussed shortly. In this queue, the service time distribution is restricted 
to be phase-type. All positive-valued distributions can be approximated by 
phase-type distributions. 

Let p be a (f-dimensional nonnegative vector whose entries sum to one, v 
be a (i-dimensional positive vector, and P be a d x d sub-stochastic matrix. 
We assume that the diagonal entries of P are zero and P is transient, namely, 
I—P is invertible. Consider a continuous-time Markov chain with d+1 phases 
(or states) where phases 1, . . . ,d are transient and phase d+ 1 is absorbing. 
For j = 1, . . . , d, the Markov chain starts in phase j with probability pj. The 
amount of time it stays in phase j is exponentially distributed with mean 
When it leaves phase j, the Markov chain enters phase t = 1, . . . ,d 
with probability Pj£ or enters phase d+1 with probability 1 — Yle=i Pji' The 
phase-type distribution with parameters (p, v, P) is the distribution of time 
from starting until absorption in phase d+1 for the above Markov chain. In 
particular, when P is a zero matrix, the associated phase- type distribution 
is a hyperexponential distribution with d phases. 

In a GI / Ph/n+GI queue, there are n identical servers working in parallel. 
The customer arrival process is a renewal process. Upon arrival, a customer 
enters service immediately if an idle server is available. Otherwise, he waits 
in a buffer with infinite waiting room that holds a first-in-first-out (FIFO) 
queue. The service times form a sequence of iid random variables, following 
a phase-type distribution. When a server finishes serving a customer, the 
server takes the leading customer from the waiting buffer. When the buffer 
is empty, the server begins to idle. Each customer has a patience time. The 



NUMERICAL ANALYSIS OF DIFFUSION MODELS 



15 



patience times are iid following a general distribution. When a customer's 
waiting time in queue exceeds his patience time, the customer abandons the 
system with no service. 

Let A be the arrival rate and 1/fi be the mean service time. The system 
is assumed to operate in the QED regime, i.e., both the arrival rate A and 
the number of servers n are large, while the traffic intensity p = A/(n/x) is 
close to one. Because customer abandonment is allowed, it is not necessary 
to assume p < 1 for the system to reach a steady state. For future purposes, 
we put 

(4.1) /3 = ^(1-P). 

Assume that the phase-type service time distribution has parameters 
(p,i/,P). Each service time can be decomposed into a number of phases. 
When a customer is in service, he must be in one of the d phases. Let Zj(t) 
denote the number of customers in phase j service at time t. In steady-state, 
one expects that the customers in service are distributed among the d phases 
following a distribution 7, given by 

(4.2) ~/ = pR- 1 p and R = (I - P') diag(i/). 

One can check that X^j=i 7j = 1 an d lj is interpreted to be the fraction of 
phase j service load on the n servers. 

Suppose that all customers, including those initial customers waiting in 
the buffer at time zero, sample their first service phases following distribution 
p upon arrival. One can stratify customers in the waiting buffer according 
to their first service phases. For j = l,...,d, we use Wj{t) to denote the 
number of waiting customers at time t whose service begins with phase j. 
Then, 

(4.3) Y j (t) = Z j (t) + W j (t) 

is the number of phase j customers in the system, either waiting or in service. 
Let Y(t) be the corresponding d-dimensional random vector and 

(4.4) Y{t) = ±={Y{t)- ni ). 

In each diffusion model, the process Y = {Y(t) : t > 0} is approximated by 
a (i-dimensional diffusion process. 
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4.2. System equation. The Gl/Ph/n + GI queue is driven by several 
primitive processes. Let E = {E(t) : t > 0} be the arrival process, where 
E(t) is the number of arrivals by time t. For j = 1, . . . , d, let Sj = {Sj(t) : 
t > 0} be a Poisson process with rate Uj, and <j)j = {4>j(i) ■ i £ N} be 
a sequence of iid (i-dimensional random vectors such that <pj(i) takes e~ 
with probability Pj£ and takes a zero vector with probability 1 — Yle=i Pjt- 
Similarly, let 4>q = {<t>a(i) '■ i £ N} be a sequence of iid d-dimensional random 
vectors such that takes e l with probability p£. For j = 0, . . . , d, define 
the routing process $j = {$j(k) : k £ N} by 

*i(*!)=x;^( 1 ')- 

2=1 

We assume that 1^(0), E, Si, ■ ■ . , S^, <&o> ■ ■ • , &d are mutually independent. 

For j = l,...,d, let Tj(t) be the cumulative amount of service effort 
received by customers in phase j service by time t. Clearly, 

(4.5) Tj(t)= [ Zj(s)ds fori>0. 

Jo 

Thus, Sj(Tj(t)) is equal in distribution to the cumulative number of phase j 
service completions by time t. (For more details, please refer to Section 4.1 
of Dai, He and Tezcan (2010) on a perturbed system.) Let Lj(t) be the 
cumulative number of phase j customers who have abandoned the system 
by time t, and L(t) be the corresponding d-dimensional vector. One can 
check that the process Y = {Y(t) : t > 0} satisfies the following equation 

d 

(4.6) Y(t) = Y(0) + * (E(t)) + E *i( S i( T i®)) - S(T(t)) - L(t), 

3=1 

where S(T(t)) = (5i(Ti(*)), . . . , S d (T d (t)))'. 

To derive the diffusion models, consider a scaled version of (4.6). We define 
several scaled processes by 

E{t) = -i=(E(t) - Xt), S(t) = -}=(S(nt) - nut), 



n yjn 

Z(t) = ^=(Z(t)-n 7 , L{t) = ^=L(t),) 
i/n Jn 



ynt\ ynt\ 



«=i v i=i 
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for t > and j = 1, . . . , d, where p 3 is the jth. column of P' . By (4.1)-(4.5), 
the dynamical equation in (4.6) turns out to be 

(4.7) 

Y(t) = y (0) - fopt + P E(t) + $o (^) 

+ £ *j (^|^) - (/ - ^)s(^) -a[' z(») i> - £(«)• 

3=1 1/0 

4.3. Diffusion models. In both diffusion models, we replace the scaled 
primitive processes in (4.7) by certain Brownian motions. These approxi- 
mations can be justified by the functional central limit theorem. When the 
number of servers n is large, the corresponding diffusion process in each 
model can be proved close to Y via a continuous map. Please refer to Dai, 
He and Tezcan (2010) for related convergence results. 

Let Be be a. one-dimensional driftless Brownian motion with variance 
Ac^/n, where (? a is the squared coefficient of variation for the interarrival 
time distribution. Let Bq, . . . , P>d, and B$ are d-dimensional driftless Brow- 
nian motions with covariance matrices H°, . . . , H d , and diag(z^), respectively, 
where 

T o _ \vk0- ~Pi) ifk = £, 
—PkPe otherwise 



and 



Hj = \P jk (l-P j£ ) ifk = £, 
~PjkPji otherwise 



for j = 1, . . . , d. We assume that Y(0), Be, Bq, . . . , Bd, Bg are mutually 
independent. In the diffusion models, the above Brownian motions take the 
places of the scaled primitive processes E, $o 5 • • • > S, respectively. Let 
Q(t) be the queue length or the number of waiting customers at time t and 

Q(t) = -^Q(t). 
Jn 



Then, Q(t) = (e'Y(t) — n) + or equivalently, 
(4.8) Q(t) = (e'Y(t)) + . 

When n is large, these waiting customers are approximately distributed 
among the d phases according to distribution p (see Lemma 2 of Dai, He 
and Tezcan (2010)), i.e., 

Wj(t) ~ Pj Q(t) forj = l,...,d. 
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It follows from (4.3) that 

Z(t)^Y(t)-pQ(t). 
By (4.4) and (4.8), this approximation has a scaled version 
(4.9) Z(t)^Y(t)-p(e'Y(t)) + . 

Let G(t) be the cumulative number of abandoned customers by time t and 

G(t) = -±=G(t). 
\ n 



These abandoned customers are also approximately distributed among the 
d phases by distribution p, i.e., 

(4.10) L(t)ttpG{t). 

We also exploit the following approximations 
(4.11) 

E(t) Xt T(t) . 1N SAT At)) . „ 

n n n n 

The approximations in (4.9)-(4.11) are used in both diffusion models. These 
two models differ only in how to approximate the scaled abandonment pro- 
cess G = {G(t) :t>0}. 

4.3.1. Diffusion model using the patience time density at zero. In the 
first diffusion model, the patience time distribution is used only through its 
density at zero when approximating G. Let F be the distribution function 
of the patience times. We assume that 

F(t) 

(4.12) F(0) = and a = lim < oo. 

Thus, a is the patience time density at zero. In this model, G is approximated 
by 

(4.13) G(t) rj a [ (e / y(s))+ds for t > 0. 

J o 

When a = 0, this approximation yields G(t) ~ for all t > 0. In this case, 
the diffusion model is for a Gl/Ph/n queue without abandonment. When 
a > 0, the intuition of (4.13) is as follows. For a many-server queue in the 
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QED regime, the queue length is typically in the order of Ofa 1 / 2 ). As a 
result, each customer's waiting time should be in the order of 0{n~ l l 2 ), 
which is relatively short because n is large. At time s, a waiting customer 
will abandon the system during the next d time units with probability a5. 
Hence, the instantaneous abandonment rate of the system is approximately 
aQ(s). It follows that 

G(t) ^ a f Q(s)ds, 
Jo 

which, along with (4.4) and (4.8), leads to the approximation in (4.13). This 
relationship can be justified by Theorem 1 of Dai and He (2010). As pointed 
out by Dai and He (2010) and Mandelbaum and Momcilovic (2009), the 
performance of a many-server queue in the QED regime is insensitive to 
the patience time distribution as long as its density a at zero is fixed and 
positive. 

In the dynamical equation (4.7), when the scaled primitive processes are 
replaced by appropriate Brownian motions and the approximations in (4.9)- 
(4.13) are employed, we obtain the following stochastic differential equation 

d 

X(0) - Pwt + P B E (t) + B (j>f£t) + B i(iP A 1)K,-7j*) 

- (J - P')B s {{p A 1)7*) - R [ (X(s) - p(e'X{s)) + ) ds 

Jo 

-pa [ (e'X(s)) + ds 
Jo 

where we take X(0) = 1^(0). We may write (4.14) into the standard form 

X(t)=X(0)+ I b{X{s))ds+ I a{X{s))dB(s), 
Jo Jo 

where for each x G the drift coefficient b is 
(4.15) b(x) = -/3/j,p - R(x - p{e'x) + ) - pa(e'x) + , 

the diffusion coefficient a is a d x d constant matrix satisfying 
a(x)a'(x) 

+ (p A 1) ( ^j Hj + ( J ~ P ') diag(^) diag( 7 )(/ - P)\ , 

V 3=1 J 



X(t) = 

(4.14) 



S(x) = 

(4.16) 
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and B is a d-dimensional standard Brownian motion. One can check that 
is positive definite and thus satisfies (2.3). The drift coefficient b in 
(4.15) is a piecewise linear function of x. Both b and a are Lipschitz con- 
tinuous. Therefore, a strong solution to (4.14) exists and is known as a 
d-dimensional piecewise OU process. In this model, the diffusion process X 
depends on the patience time distribution only through its density at zero. 
When using the proposed algorithm to solve the stationary density, it follows 
from Proposition 2 that the linear system (3.12) has a unique solution. 

If we replace p by one in (4.14), the resulting diffusion process turns out 
to be the diffusion limit for G/Ph/n + GI queues in Theorem 2 of Dai, He 
and Tezcan (2010). This limit process is the strong solution of the following 
stochastic differential equation 



Since p is close to one in the current setting, the above limit process justifies 
the diffusion model in (4.14). 

4.3.2. Diffusion model using patience time hazard rate scaling. When 
the patience time distribution does not have a density at zero, the diffusion 
model in (4.14) fails to exist. When a = and p > 1, the diffusion process 
X in (4.14) does not have a stationary distribution. In this case, the model 
cannot be a satisfactory approximation of the many-server queue, as the 
queue may have a stationary distribution thanks to customer abandonment. 
It is also demonstrated in Reed and Tezcan (2009) that when the density near 
zero rapidly changes, the system performance can be sensitive to the patience 
time distribution in a neighborhood of the origin. In this case, using the 
patience time density at zero solely may not yield adequate approximation 
to the queue. Our second diffusion model exploits the idea of scaling the 
patience time hazard rate function, which was first proposed in Reed and 
Ward (2008) for single-server queues and was recently extended to many- 
server queues in Reed and Tezcan (2009). 

In this model, we assume that the patience time distribution F satisfies 



d 



X{t) = X(0) - (3ppt+pB E (t) + B (pt) +^Bj{v nj t) 



(4.17) 




F(0) = 
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and it has a bounded hazard rate function h, given by 

hit) = for t > 0, 

where fp is the density of F. With the hazard rate function, F can be 
written by 

F(t) = 1 - exp ( - J h(s) da) for t > 0. 

In the second diffusion model, the scaled abandonment process G is approx- 
imated by 

(4.18) G(t)« J J h(^—)duds fori>0. 

The entire patience time distribution is built into the approximation through 
its hazard rate function. The intuition of the hazard rate scaling approxima- 
tion was explained in Reed and Ward (2008): Consider the Q(s) waiting cus- 
tomers in the buffer at time s. In general, only a small fraction of customers 
can abandon the system when the queue is working in the QED regime. 
Then by time s, the ith customer from the back of the queue has been 
waiting around i/X time units. Approximately, this customer will abandon 
the queue during the next 5 time units with probability h(i/X)8. It follows 
that for the system, the instantaneous abandonment rate at time s is close 
to Ylf=x By (4.4) and (4.8), the scaled abandonment rate can be 

approximated by 

^X>(aW HV) d » = y "(V) d "' 

from which (4.18) follows. Note that the arrival rate A is in the order of 
0{n) and Q(s) is in the order of 0(n 1//2 ). The patience time distribution 
in a small neighborhood of zero, not just its density at zero, is considered 
in the instantaneous abandonment rate in (4.19). Hence, the hazard rate 
scaling approximation in (4.18) is more accurate than that in (4.13). This 
approximation can be justified for GI/M/n + GI queues by Propositions 9.1 
and 9.2 in Reed and Tezcan (2009). With minor modifications to the proofs, 
these two propositions can be extended to Gl/Ph/n + GI queues. 

Let m be a nonnegative integer. Suppose that the hazard rate function h 
is m times continuously differentiable in a neighborhood of zero. By Taylor's 
theorem, 



h{z) « fc(0) + jTfrM(O) 
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for z > small enough, where is the ith. order derivative of h. In this 
case, the approximation in (4.18) turns out to be 

G(t) » h(0) j\e'Y(s)) + ds + £j ^ffiff j\(e'Y (s)) + Y +1 ds. 

Because /i(0) is identical to the patience time density at zero, the approxima- 
tion in (4.13) can be regarded as the zeroth degree Taylor's approximation 
of (4.18). When the patience times are exponentially distributed, the hazard 
rate function is constant and the two approximations in (4.13) and (4.18) 
are identical. See Section 4 of Reed and Tezcan (2009) for more discussion. 

Using the Brownian motion replacement and the approximations in (4.9)- 
(4.11) and (4.18), we obtain the second diffusion model for the Gl/Ph/n + 
GI queue, given by 



X(t) = X(0) - f3fipt + P B E (t) + Bo{pnt) + Y, B A(P A 1)^77*) 
( 4 - 2 °) - (/ - P')B s ((p A 1)t*) - R J\x{s) - p(e'X(s))+) ds 



J 



The diffusion process X in (4.20) has the same diffusion coefficient a as in 

the first model (4.14). Its drift coefficient b is 

(4.21) 

r(e'x)+ */nu\ j 
b( x ) = -Ppp- R(x -p{e'x) + ) -p j hl^^-jdu for x G R . 

Because h is bounded, the drift coefficient b is Lipschitz continuous and 
the stochastic differential equation (4.20) has a strong solution. By Propo- 
sition 2, the solution to the linear system (3.12) is unique when we use the 
proposed algorithm to solve the stationary density of this diffusion model. 
Comparing (4.14) and (4.20), one can see that the two models differ only in 
how the patience time distribution is incorporated. Because a more accurate 
approximation is used for the abandonment process, the second model can 
provide a better approximation for the queue. 



5. Choosing a reference density. The reference density controls the 
convergence of the proposed algorithm. In this section, we discuss how to 
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choose appropriate reference densities for the diffusion models. Some con- 
siderations are as follows. 

First, to be a reference density, a candidate function r must satisfy (3.1) 
even though the stationary density g is unknown. The second condition in 
(3.1) requires that r have a comparable or slower decay rate than g. When g 
is bounded, its decay rate is sufficient to determine a function r that satisfies 
(3.1). 

Second, the most computational effort in the algorithm is constructing and 
solving the system of linear equations (3.12). As demonstrated by Proposi- 
tion 3, the finite-dimensional approximates the infinite-dimensional H 
better as k increases, thus reducing the approximation error. On the other 
hand, as the dimension of Hk increases, constructing and solving (3.12) re- 
quires more computation time and memory space. The condition number 
of the matrix A in (3.12) also gets worse as the dimension of becomes 
large. This yields higher round-off error. A "good" reference density should 
balance the approximation error and the round-off error. With such a ref- 
erence density, it is possible to have small approximation error even if the 
dimension of is moderate. 

Intuitively, when r is "close" to the stationary density g, both the ratio 
function q and the projection c are close to constant. We can thus expect 
that a space with a moderate dimension is able to produce a satisfac- 
tory approximation. All these observations motivate us to explore the tail 
behavior of a diffusion model. 

5.1. Tail behavior. Let us focus on the diffusion limit in (4.17). Assume 
that the piecewise OU process X is positive recurrent and has a stationary 
distribution. Let X{po) be the corresponding <i-dimensional random vector 
in steady state. Since p is close to one, the tail behavior of the diffusion 
process X in (4.14) is expected to be comparable to that of the limit diffusion 
process X in (4.17). 

To explore the tail behavior of X(oo), consider a sequence of GI/GI/n + 
GI queues in the QED regime. If all patience times are infinite, the queues 
turn out to be GI/GI/n queues without customer abandonment. In each 
queue, the service times are iid following a general distribution. We assume 
that these queues, each indexed by the number of servers n, have the same 
service time distribution. Let A n be the arrival rate of the nth system. To 
mathematically define the QED regime, we assume that 



(5.1) 
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and 



(5.2) 



lim \/n(l — p n ) = ft for some $ 6 

n— >oo 



where p n = \ n / (nu) is the traffic intensity of the nth system. 

Assume that all these queues are in steady state. Let N n (oo) be the sta- 
tionary number of customers in the nth system and 



iVn(oc) 



1 



n 



(N n (oo)-n). 



For GI/GI/n queues in the QED regime, the limit queue length in steady 
state was studied in Gamarnik and Momcilovic (2008), where the service 
time distribution is assumed to be lattice- valued on a finite support. The 
authors first showed that N n (oo) weakly converges to a random variable 
N(oo) as n goes to infinity, and then proved that 



(5.3) 



lim 1 logP[A>(oo) 



> Z 



2 i 2 ' 



where (? a and (? s are the squared coefficients of variation of the interarrival 
and the service time distributions, respectively. In (5.3), the decay rate does 
not depend on the service time distribution beyond its first two moments. 
Recently, this result has been extended in Gamarnik and Goldberg (2011) 
to GI/GI/n queues with a general service time distribution. 

When a = and d = 1, the limit diffusion process X in (4.17) is for 
GI/M/n queues without customer abandonment. In this case, the service 
time distribution is exponential and N(oo) = X(oo). It was proved in Halfin 
and Whitt (1981) that the stationary density of X(oo) has a closed-form 
expression 



(5.4) 



9W 



a± exp 
ci2 exp 



2/3z 



if z < 0, 
if z > 0, 



where a\ and a>2 are normalizing constants making g continuous at zero. The 
decay rate of g in (5.4) is consistent with (5.3). Both formulas suggest that 
iV(oo) has an exponential tail on the right side. 

For a GI/GI/n + GI queue with many servers and customer abandon- 
ment, the limiting tail behavior of N n (oo) remains unknown except for very 
simple cases. When a > and d = 1, the limit diffusion process X in (4.17) 
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is a one-dimensional piecewise OU process. It admits a piecewise normal 
stationary density 



(5.5) g{z) 



ct 3 exp ( i + c 2 I if 2 < 0, 



where a% and 04 are normalizing constants that make g continuous at zero. 
See Browne and Whitt (1995) for more details. 

Observing (5.3) and (5.5), we conjecture that for a sequence of GI/GI/n+ 
GI queues in the QED regime, the limiting tail behavior of N n (oo) depends 
on the service time distribution only through its first two moments, and on 
the patience time distribution only through its density at zero. 

Conjecture 3. Consider a sequence of GI/GI/n + GI queues that 
satisfies (4-12), (5.1), and (5.2). Assume that the patience time distribution 
has a positive density at zero, i.e., a > in (4. 12). Assume further that 
the interarrival and the service time distributions satisfy the Tq assumptions 
(i)-(iii) in Section 2.1 of Gamarnik and Goldberg (2011). Then, (a) N n (oo) 
exists for each n; (b) the sequence of random variables {N n (oo) : n G N} 
weakly converges to a random variable iV(oo); (c) iV(oo) satisfies 

1 - a 

lim ^logP[iV(oo) > z] 



z^oo z 2 " ' /i(cg + cg)' 

The intuition below may help understand why the conjectured decay rate 
must be Gaussian. When N(oo) > z for some z > 0, there are more than 
n l / 2 z waiting customers in the queue correspondingly, and each waiting cus- 
tomer is "racing" to abandon the system. At any time, the instantaneous 
abandonment rate is approximately proportional to the queue length. In 
such a system, the customer departure process, including both service com- 
pletions and customer abandonments, behaves as if the system is a queue 
with infinite servers. Thus, one can expect that the tail of the limit queue 
length is Gaussian, which decays much faster than an exponential tail for 
queues without abandonment. 

5.2. Reference densities for model (4-14)- For Gl/Ph/n + GI queues, 
the limit diffusion process X in (4.17) satisfies 

iV(oo) = e'X{oo). 

The discussion in Section 5.1 gives ample evidence of the tail behavior 
P[iV(oo) > z] as z -> 00. Although the left tail P[iV(oo) < -z] as z -> 00 
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remains unknown when d > 1, our numerical experiments suggest that this 
tail is not sensitive to the service time distribution beyond its mean. Thus, 
we use the left tail for a queue with an exponential service time distribution 
to construct the reference density. We propose to use a product reference 
density 

d 

(5.6) r(x) = Yl rj{xj) for x G R d . 

j=l 

When a = and p < 1 in (4.14), there is no abandonment in the queue. 
Based on (5.3) and (5.4), we choose 



(5.7) r 3 {z) 



eX P I " .2 " ) lf Z ^ °> 



d + cj i + c 



where f3 is given by (4.1). The function rj has an exponential tail on the right 
and a Gaussian tail on the left. One can check that the reference density 
given by (5.6) and (5.7) satisfies condition (3.3). In (5.7), we set the shift 
term for z < to be 7j/3 according to the following observation. In the 
associated queue, /? is the scaled mean number of idle servers and jj is the 
fraction of phase j service load. In steady state, one can expect that Yj(t), 
the centered and scaled number of phase j customers, is around —Jjfi. 

When a > in (4.14), the associated queue has abandonment. By (5.5) 
and Conjecture 3, we choose 
(5.8) 

exp(-^) 

whose two tails are both Gaussian but have different decay rates. This refer- 
ence density also satisfies (3.3). In (5.8), the shift term for z > is taken to 
be pj/i/3/a because of the observation below. When p > 1, the throughput of 
the queue is nearly n\x. Let go be the scaled queue length "in equilibrium", 
i.e., the arrival and the departure rates of the system are balanced when the 
queue length is around n 1 / 2 ^- Because in this case the abandonment rate is 
an 1 / 2 ^, we must have A = np + an l / 2 qQ, or qo = —pf3/a by (4.1). Since the 
fraction of phase j waiting customers is around pj, Yj(t) is around —pjp,f3/a 
as the queue reaches a steady state. 
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5.3. Reference densities for model (4-20). For the diffusion model (4.20) 
that adopts the patience time hazard rate scaling, the tail behavior of X 
in steady state is left to future research. In some cases, we may exploit the 
diffusion limit in (4.17) to facilitate the choice of a reference density for the 
current model. The principle is again to ensure that the reference density 
has a comparable or slower decay rate than the stationary density of X. For 
that, we build an auxiliary queue that shares the same arrival process and 
service times with the Gl/Ph/n + GI queue, but the auxiliary queue may 
have no abandonment or have an exponential patience time distribution. 
Let X be the diffusion process in (4.14) for the auxiliary queue. If X has 
a slower decay rate than X, a reference density of X must be a reference 
density of X, too. 

When p < 1, the auxiliary queue is a Gl/Ph/n queue. It is intuitive that 
the queue length decays faster in the Gl/Ph/n + GI queue than in the 
auxiliary queue because the latter has no abandonment. As a consequence, 
X has a slower decay rate than X and the reference density given by (5.6) 
and (5.7) for X can be used for the current model. 

When p > 1, the auxiliary queue is a Gl/Ph/n + M queue. Let a > 
be the rate of the exponential patience time distribution, which is to be 
determined in order for X to have an appropriate decay rate. For that, we 
need investigate the abandonment process of the Gl/Ph/n + GI queue. 

Assume that the hazard rate function h is m times continuously differen- 
tiable in a neighborhood of zero for some nonnegative integer m, and among 
I = 0, ...,m, there is at least one h^(0) ^ 0. We follow the convention 
that /i(°)(0) = h(0). Let £ be the smallest nonnegative integer such that 
hV°\0) ^ 0. For z > in a small neighborhood of zero, the loth, degree 
Taylor's approximation of h is 

, , ^°)(0)^° 
(5.9) h(z) « jj—, 

which, along with (4.8) and (4.18), implies that the scaled abandonment 
process can be approximated by 

G(t)K WT§ i d, 

This approximation implies that the abandonment process depends on the 
hazard rate function primarily through /i^°)(0), the nonzero derivative at 
the origin with the lowest order. It also implies that the scaled abandonment 
rate at time t is approximately 

$(*) u (s/nu\ u _ ™ V2/iW)) (0),/v^o+i 



(5.10) / MVJ d - A^o + i)! (9(t)) 
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In the hazard rate scaling, the scaled queue length in equilibrium qo sat- 
isfies 



qo > - >nu 



(5.11) A = n/i + V^/ M^~~ ) dn - 

If (5.10) holds, it turns out to be 



which gives us 
(5.12) q 



n (.to+i)/2 h (M(p) . 



1 /A'-i',, • l;!jA »;/)\ : 



V h( e o)(0) 



The scaled queue length process fluctuates around this equilibrium length. 
Correspondingly, the instantaneous abandonment rate changes around an 
equilibrium level, too. This observation motivates us to take 

(5 13) a - "^^(ty 

for the auxiliary GI / Ph/n + M queue. With this setting, the original queue 
and the auxiliary queue have comparable abandonment rates when the 
scaled queue length is close to qo. For any qi > qo, when the scaled queue 
length is q\ in both queues, the abandonment rate in the auxiliary queue is 
lower because 

n to/2 h (M(Q) . 

aqi < A*(/o + l)ffi • 
Hence, when the queue length is longer than qo, it decays slower in the 
auxiliary queue than in the original queue. Consequently, the decay rate of 
X is slower than that of X and the reference density of X can work for this 
diffusion model. 

The above discussion suggests a product reference density in (5.6) with 



(5.14) Tj{z) 



exp ( - - ; ^ 2 ) if z < 0, 



exp( _ °(Z-PJ9>) 2 + <*M _ -jg., ifz>Q 



where qo follows (5.12) and a follows (5.13). 

The above reference density fails when p = 1 and lo > 0, because qo = 
by (5.12) and thus a is zero in (5.13). In this case, we can still choose a 
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reference density by (5.6) and (5.14) but using a traffic intensity p that is 
slightly larger than one. Because the tail of the queue length becomes heavier 
as p increases, a reference density for model (4.20) with p > 1 must have a 
comparable or slower decay rate than the stationary density of the model 
with p = 1. 

The reference density in (5.6) and (5.14) that exploits the lowest-order 
nonzero derivative at the origin may fail when the hazard rate function has 
a rapid change near the origin. In this case, the Taylor's approximation in 
(5.9) may not be satisfactory when the queue length is not short enough. 
Such an example is discussed in Section 6.4. Choosing a reference density 
for that is more flexible. In addition, the above procedure cannot choose a 
reference density when all /i^(0)'s are zero, i.e., the hazard rate function is 
zero in a neighborhood of the origin. This topic will be explored in future 
research. 

5.4. Truncation hypercube. Once the reference density has been deter- 
mined, we can choose the truncation hypercube K in (3.14) by the procedure 
below. First, pick a small number eo > 0. Then, choose a hypercube K such 
that 



JR d \K 

When Eo is small enough, the influence of the reference density outside K is 
negligible in computing the stationary density. 

6. Numerical examples. Several numerical examples are presented in 
this section. In each example, we compute the stationary distribution of the 
number of customers in a many-server queue using a diffusion model and 
the proposed algorithm. We assume that the customer arrivals follow a ho- 
mogeneous Poisson process and the service times follow a two-phase hyper- 
exponential distribution with mean one, i.e., the system is an M / H2/n + GI 
queue with c\ = 1 and p = 1. In such a queue, there are two types of 
customers. One type has a shorter mean service time than the other, and 
the service times of either type are iid following an exponential distribu- 
tion. We approximate this queue by a two-dimensional diffusion process X. 
When the patience time distribution is exponential, both (4.14) and (4.20) 
yield the same diffusion process. When the patience time distribution is 
non-exponential, we use model (4.20) as it is more accurate. The results 
computed using the diffusion models are compared with the ones obtained 
either by the matrix-analytic method or by simulation. Please refer to Neuts 
(1981) and Latouche and Ramaswami (1999) for the implementation of the 



(5.15) 
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-diffusion model I 
-matrix-analytic 




the number of customers in system 

(a) p = 1.141 and n = 50 




the number of customers in system 

(b) p = 1.045 and n = 500 



Fig 1: The stationary distribution of the customer number in the M / H2/n + 
M queue. 



matrix-analytic method. All simulation results are obtained by averaging 
twenty runs and in each run, the queue is simulated for one hundred thou- 
sand time units. 

In the proposed algorithm, all numerical integration is implemented using 
a Gauss-Legendre quadrature rule. See Kress (1998) for more details. When 
computing An or Vi in (3.13), the integrand is evaluated at eight points in 
each dimension. In the numerical examples, the tail probability 

(6.1) P[A"i(oo)+X 2 (oo) > z] = / g{x) dx for some z£R 

J {x£R 2 :xi+X2>z} 

is also computed, where X(oo) = (Xx(oo), ^(oo))' is the two-dimensional 
random vector having probability density g. The integral in (6.1) is com- 
puted by adding up the integrals over the finite elements that intersect with 
the set {x € M 2 : x\ + X2 > z}, and the integral over each finite element 
is again computed using a Gaussian-Legendre quadrature formula. Because 
the indicator function has jumps inside certain finite elements, we use sixty- 
four points in each dimension when evaluating the integrand over each finite 
element. 

6.1. Example 1: an M/H2/T1 + M queue. Consider an M/H2/n + M 
queue that has an exponential patience time distribution. We are interested 
in such a queue because its customer-count process iV = {N(t) : t > 0} 
is a quasi-birth-death process, where N(t) is the number of customers in 
system at time t. The stationary distribution of that can be computed by 
the matrix-analytic method. 
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Model (4.14) 


Matrix-analytic 


Mean queue length 


17.27 


17.16 


Abandonment fraction 


0.1512 


0.1503 


P[JV(oo) > 45] 


0.8675 


0.8523 


P[iV(oo) > 50] 


0.6785 


0.6726 


P[iV(oo) > 100] 


0.08700 


0.07436 


P[AT(oo) > 130] 


0.008662 


0.003299 


(a) p = 1.141 and n — 50 






Model (4.14) 


Matrix-analytic 


Mean queue length 


54.17 


54.05 


Abandonment fraction 


0.05181 


0.05173 


P[iV(oo) > 470] 


0.9701 


0.9694 


P[AT(oo) > 500] 


0.6838 


0.6818 


P[iV(oo) > 600] 


0.2244 


0.2229 


P[iV(oo) > 750] 


0.008233 


0.006395 



(b) p = 1.045 and n = 500 
Table 1 

Performance measures of the M/H2/11 + M queue. 



In this example, we take a = 0.5 for the rate of the exponential patience 
time distribution and 

p = (0.9351, 0.0649)' and v = (9.354, 0.072)' 

for the hyperexponential service time distribution. The mean service time of 
the second-type customers is more than one hundred times longer than that 
of the first type. Although over ninety percent of customers are of the first 
type, the fraction of its workload is merely ten percent, i.e., 7 = (0.1,0.9)'. 
Such a distribution has a large squared coefficient of variation c 2 s = 24. 

The queue is approximated by the two-dimensional piecewise OU process 
X in (4.14). Because the service time distribution is hyperexponential, P is a 
zero matrix and thus R = diag(z^). By (4.15) and (4.16), the drift coefficient 
of X is 

(6 2) b(x) = (~ PlfJ, P ~ Vl ^ Xl ~P^ Xl + x ?) + ) -P\a(xi +x 2 ) 
' ' '"' - J/ 2(x2 ~P2{x\ +x 2 ) + ) -p 2 a(xi +x 2 ) 

and the covariance matrix of the diffusion coefficient is 

W(p+(Ml)) 



(6 - 3) 1 A 1)) 



for all x £R 2 . 
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0.08 



20 40 60 

the number of customers in system 



80 



Fig 2: The stationary distribution of the customer number in the M / H2/n + 
M queue, with p = 1.112 and n = 20. 

Three scenarios are considered in this example, in all of which the queue 
is overloaded. In the first two scenarios, there are n = 50 and 500 servers, 
respectively. The arrival rates are A = 57.071 and 522.36, or equivalently, 
p = 1.141 and 1.045. By (4.1), (3 = —1 in both scenarios. The third scenario, 
with n = 20 servers, will be presented shortly. 

To compute the stationary distribution of X, we use a product reference 
density given by (5.6) and (5.8). To generate basis functions by the finite 
element method, we set the truncation rectangle K = [—7,32] x [—7,32], 
which is obtained by (5.15) with eq = 10 -7 , and use a lattice mesh in which 
all finite elements are 0.5 x 0.5 squares. 

Once the stationary density of X is obtained, one can approximately 
produce the distribution of iV(oo), the stationary number of customers in 
system. Note that the probability density of X\{po) + A^oo) is given by 



For the first two scenarios, the distributions of N(oo) obtained by the dif- 
fusion model are illustrated in Figure 1. In the same figure, the stationary 
distributions computed by the matrix-analytic method are plotted, too. We 
see good agreement in Figure 1. Comparing the two scenarios, we also find 
out that the diffusion model in (4.14) is more accurate when the number of 




for zeR. 



The distribution of iV(oo) can be approximated by 




1 



for % = 0, 1, . . . 
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servers n is larger. This observation is consistent with the many-server limit 
theorem for G/Ph/n + GI queues in Dai, He and Tezcan (2010). 

The matrix- analytic method can be used because in this queue, the three- 
dimensional process {(Q(t), Z\(t), Zzit)) : t > 0} forms a continuous-time 
Markov chain and the customer-count process N is a quasi-birth-death pro- 
cess. Clearly, N(t) = Q(t) + Z\(t) + Zi($)- At time t, N is said to be at level 
£ if N(t) = I. In this example, level £ consists of I + 1 states if £ < n and it 
contains n + 1 states if £ > n. In the matrix-analytic method, the transition 
rate matrices between adjacent levels are exploited to compute the station- 
ary distribution of N iteratively. Each iteration requires 0(n 3 ) arithmetic 
operations. For this queue, the transition rate matrices at different levels are 
different because the abandonment rate depends on the queue length. For 
implementation purposes, we assume in the algorithm that at level £ > £q 
for some £q 2> n, the abandonment rate at level £ is a(lo — n) rather than 
a(£ — n). In other words, the transition rate matrices at level I are invariant 
with respect to £ when £ > £q. We take £q = n + 2000 in all numerical 
examples. The extra error caused by this modification is negligible, because 
in this queue, the queue length is in the order of 0(n 1//2 ) and the chance of 
the customer number exceeding £$ is extremely rare. 

To investigate the diffusion model in (4.14) quantitatively, we list some 
steady-state performance measures in Table 1. They include the mean queue 
length, the fraction of abandoned customers, and the probabilities that the 
number of customers exceeds certain levels. Using the diffusion model, 

the mean queue length ~ y/n I {x\ + X2) + g{x) dx 

and 

the mean number of idle servers ~ \fn \ (x\ + X2)~ g(x) dx. 

Jr 2 

It follows from the latter approximation that 

the abandonment teaetion p, 1 - t( a - ^ f (xi + „)-,(,) ix 

A V J R 2 

In the table, the tail probability P[A r (oo) > £] is approximated by 

P[JV(oo) >£] «pfXi(oo) + X 2 (oo) > —(£-n)] fori = 0,1, 
L Jn 



and P[Xi(oo) + A^oo) > (£ — n)/y/n\ is computed via (6.1). In both sce- 
narios, the diffusion model produces satisfactory numerical estimates. 
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20 40 60 80 100 400 500 600 700 

the number of customers in system the number of customers in system 



(a) p = 0.8586 and n = 50 (b) p = 0.9553 and n = 500 

Fig 3: The stationary distribution of the customer number in the MjH^jn 
queue. 

The computational complexity of the proposed algorithm, whether in 
computation time or in memory space, does not change with the number 
of servers n. In contrast, the matrix-analytic method becomes computation- 
ally expensive when n is large. In particular, the memory usage becomes 
a serious constraint when a huge number of iterations are required. For 
the n = 500 scenario in this example, it took around one hour to finish 
the matrix-analytic computation and the peak memory usage is nearly five 
gigabytes. Using the diffusion model and the proposed algorithm, it took 
less than one minute and the peak memory usage is less than two hundred 
megabytes on the same computer. See Section 7.4 for more discussion on 
the computational complexity. 

Although the diffusion model is motivated and derived from the theory of 
many-server queues, it is still relevant for a queue with a modest number of 
servers. In the third scenario, there are n = 20 servers and the arrival rate 
is A = 22.24. Thus, p = 1.112 and /3 = —0.5. In the proposed algorithm, we 
keep the same truncation rectangle and lattice mesh as in the previous two 
scenarios, and the reference density is again from (5.6) and (5.8). As illus- 
trated in Figure 2, the diffusion model can still capture the exact stationary 
distribution for a queue with as few as twenty servers. 

6.2. Example 2: an M jH^jn queue. In this example, an M/H2/11 queue 
without abandonment is considered. The hyperexponential service time dis- 
tribution has 



p = (0.5915, 0.4085)' and v = (5.917, 0.454)'. 
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0.08 



diffusion model 

matrix- analytic 



tt 0.04 

c 



0.02 



0.06 




20 40 60 80 

the number of customers in system 



100 



Fig 4: The stationary distribution of the customer number in the MjH^jn 
queue, with p = 0.8882 and n = 20. 

Thus, c 2 s = 3 and 7 = (0.1,0.9) ; . Since there is no abandonment, we must 
take p < 1 in order for the system to reach a steady state. 

The diffusion model in (4.14) with a = is used. The drift and the 
diffusion coefficients of X are given by (6.2) and (6.3). The first scenario has 
n = 50 servers and the second scenario has n = 500 servers. The respective 
arrival rates are A = 42.929 and 477.64. Hence, p = 0.8586 and 0.9553, 
both yielding f3 = 1. The product reference density is given by (5.6) and 
(5.7). With eo = 10~ 7 , the truncation rectangle is set by (5.15) to be K = 
[—7,35] x [—7,35], which is divided into 0.5 x 0.5 finite elements. 

The stationary distribution of the number of customers in system is shown 
in Figure 3. In both scenarios, the diffusion model produces a good approx- 
imation of the result by the matrix-analytic method. As in the previous ex- 
ample, the diffusion model is more accurate when the system scale is larger. 
Several performance measures in steady state are listed in Table 2. As in 
Table 1, satisfactory agreement can be found between the two approaches. 

The third scenario has n = 20 servers with arrival rate A = 17.76. Then, 
p = 0.8882 and /3 = 0.5. With eo = 10 -7 , the truncation rectangle is taken 
to be K = [—7, 79] x [—7, 79]. The lattice mesh consists of 0.5 x 0.5 finite ele- 
ments. The distribution of N(oo) is shown in Figure 4. For a queue without 
abandonment, the diffusion model is still useful when the number of servers 
is modest. 

6.3. Example 3: an M/H2/n + E^ queue. Consider an M/H2/n + E^ 
queue, where k > 1 is a positive integer and +E^ signifies an Erlang-A; 
patience time distribution. In this queue, each patience time is the sum of k 
stages and the stages are iid having an exponential distribution with mean 
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Model (4.14) 


Matrix-analytic 


Mean queue length 


2.267 


2.419 


P[JV(oo) > 40] 


0.6908 


0.6578 


P[JV(oo) > 50] 


0.2072 


0.2012 


P[iv(oo) > 70J 


0.03395 


0.03655 


lr[iv^ooj .> 1UUJ 


u.uuooo / 




(a) p = 


0.8586 and n = 


50 




Model (4.14) 


Matrix-analytic 


Mean queue length 


8.753 


8.800 


P[A r (oo) > 450] 


0.9038 


0.9005 


P[JV(oo) > 500] 


0.2285 


0.2263 


P[iV(oo) > 600] 


0.01910 


0.01908 


P[iV(oo) > 700] 


0.002241 


0.001903 


(b)p = 


0.9553 and n — 


500 




Table 2 





Performance measures of the M/Hz/n queue. 



1/8. When k > 1, the probability density at zero of an Erlang-A; distribution 
is zero. The diffusion model in (4.14) does not have a stationary distribution 
when the queue is overloaded. Hence, we evaluate the diffusion model in 
(4.20) that exploits the patience time hazard rate scaling. In the following 
numerical experiments, we take k = 2 or 3 for the Erlang-A; distribution and 
set 6 = k, so the mean patience time is one unit time. The hyperexponential 
service time distribution is taken to be identical to that in Section 6.2. 
The hazard rate function of the Erlang-/c distribution is 

nkjk— 1 

Ht) = for t > 0. 

£=0 11 

For the diffusion model (4.20), it follows from (4.21) that the drift coefficient 
of X is 

(6 4) b(x) = (~' Px ^ ~ Vx ^ xx ~ Pl ( Xl + X2 ) + ) ~ P ir l((xi + x 2) + ) 

~ V 2 {x 2 - P2{X\ + X 2 ) + ) - P2V((X1 + ^2) + ) 



where 



'0 



,«^r^) da ^,-A log ( E ^) fora ,„. 

• m=0 



The first two scenarios has n = 50 and 500 servers, respectively. Their 
respective arrival rates are A = 42.929 and 477.64. Hence, p = 0.8586 and 
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+E 2 


+E 3 




Model (4.20) 


Simulation 


Model (4.20) 


Simulation 


Mean queue length 


0.9820 


1.061 


1.201 


1.302 


Abandonment fraction 


0.007974 


0.008592 


0.005629 


0.006115 


¥[N(oo) > 35] 


0.8881 


0.8745 


0.8896 


0.8762 


P[7V(oo) > 40] 


0.6755 


0.6399 


0.6798 


0.6448 


P[iV(oo) > 50] 


0.1671 


0.1581 


0.1788 


0.1707 


P[JV(oo) > 60] 


0.03238 


0.03353 


0.04420 


0.04584 




(a) p = 0.8586 and n = 50 








+E 2 




+E 3 






Model (4.20) 


Simulation 


Model (4.20) 


Simulation 


Mean queue length 


4.960 


5.048 


6.455 


6.569 


Abandonment fraction 


0.001689 


0.001729 


0.0007611 


0.0007931 


P[iV(oo) > 450] 


0.9003 


0.8964 


0.9022 


0.8984 


P[iV(oo) > 480] 


0.4759 


0.4643 


0.4859 


0.4746 


P[JV(oo) > 500] 


0.1995 


0.1966 


0.2151 


0.2124 


P[JV(oo) > 550] 


0.02798 


0.02841 


0.04412 


0.04458 



(b) p = 0.9553 and n = 500 
Table 3 

Performance measures of the MjHiju + Ek queue with p < 1. 



0.9553, both leading to (3 = 1. In the proposed algorithm, the reference 
density is chosen according to (5.6) and (5.7). The truncation rectangle is 
taken to be K = [—7,35] x [—7,35] and is divided into 0.5 x 0.5 finite 
elements. Some performance estimates can be found in Table 3. 

The third and fourth scenarios are for the case p > 1. They have n = 50 
and 500 servers, and arrival rates A = 57.071 and 522.36, respectively. Then, 
p = 1.141 and 1.045, both having [3 = —1. For these two scenarios, we adopt 
the reference density in (5.6) and (5.14). When k = 2, each patience time 
has two stages. The hazard rate function of the patience time distribution 
has h(0) = and h^(0) = 9 2 , so £ = 1 in (5.12) and (5.13). Because 
a in (5.13) depends on n, both the reference density and the truncation 
rectangle change with n. With eo = 10 -7 , the truncation rectangle is set 
to be K = [-7, 13] x [-7, 13] for n = 50 and to be K = [-7, 16] x [-7, 16] 
for n = 500. When k = 3, a patience time consists of three stages. In 
this case, h(0) = h^(0) = and h^{0) = 86> 3 , so £ = 2. We set K = 
[-7,11] x [-7,11] for n = 50 and K = [-7,15] x [-7,15] for n = 500. 
All truncation rectangles are partitioned into 0.5 x 0.5 finite elements. The 
performance estimates are listed in Table 4. 

To evaluate the diffusion model (4.20), we list corresponding simulation 
estimates of the performance measures in both tables. As in the previous 
examples, the diffusion model produces adequate performance approxima- 
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+E 2 


+E 3 




Model (4.20) 


Simulation 


Model (4.20) 


Simulation 


Mean queue length 


15.03 


14.94 


19.44 




19.31 


Abandonment fraction 


0.1332 


0.1334 


0.1303 




0.1305 


P[A r (oo) > 45] 


0.9568 


0.9490 


0.9704 




0.9645 


P[AT(oo) > 50] 


0.8780 


0.8648 


0.9169 




0.9066 


P[AT(oo) > 70] 


0.3325 


0.3121 


0.5037 




0.4761 


P[JV(oo) > 90] 


0.008153 


0.009354 


0.03033 




0.03422 




(a) p = 1.141 and n — 50 










+E 2 






+E 3 






Model (4.20) 


Simulation 


Model (4.20) 


Simulation 


Mean queue length 


76.50 


76.20 


119.5 




119.1 


Abandonment fraction 


0.04438 


0.04437 


0.04340 




0.04337 


P[Ar(oo) > 480] 


0.9857 


0.9846 


0.9946 




0.9940 


P[AT(oo) > 500] 


0.9390 


0.9363 


0.9770 




0.9756 


P[N{oo) > 600] 


0.3115 


0.3051 


0.6733 




0.6645 


P[AT(oo) > 700] 


0.0009757 


0.0009658 


0.04260 




0.04358 




(b) p = 1.045 and n = 500 









Table 4 

Performance measures of the Ml Hi In + Ek queue with p > 1. 

tions. 

Theoretically, the matrix-analytic method can be used in this example 
as the customer-count process N is also a quasi-birth-death process. But 
it is impractical because the computational complexity is too high. Con- 
sider the case k = 2. Let Vi(t) and V^(t) be the respective numbers of 
waiting customers whose patience times are in the first and in the second 
stage at time t. For this M/H^/n + E2 queue, the four-dimensional process 
{(Vi(t),V2(t), Z\(t), Z2(t)) : t > 0} is a continuous-time Markov chain. At 
level £, there are £ + 1 states if £ < n and there are (re + 1) (£ — n + 1) states 
if t > n. The number of states at level £ is formidable when £ is large. Even 
if we may truncate the state space using the technique described in Sec- 
tion 6.1, the number of states is still too large to apply the matrix-analytic 
method. In fact, we are not aware of any other numerical methods other 
than simulation that can produce approximations in Tables 3 and 4. 

6.4. Example 4-' o- n MlH^jn + H2 queue. Let us consider an example 
in which the patience time hazard rate function changes rapidly near the 
origin. As pointed out by Reed and Tezcan (2009), the performance of such 
a queue is sensitive to the patience time distribution in a neighborhood of 
zero. A model that exploits the patience time density at zero solely may 
not produce adequate performance estimates. In this example, the patience 
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Model (4.14) 


Model (4.20) 


Simulation 


Mean queue length 


0.4709 


4.869 


4.845 


Abandonment fraction 


0.1714 


0.1504 


0.1499 


P[JV(oo) > 40] 


0.9578 


0.9749 


0.9728 


P[iV(oo) > 50] 


0.3158 


0.6377 


0.6111 


P[A(oo) > 60] 


1.044 x 10~ 7 


0.1895 


0.1737 


P[A(oo) > 70] 


1.097 x 10 


0.02568 


{ \ a A 1 A C\ 

0.02142 


(a 


p = 1.141 and n 


= 50 






Model (4.14) 


Model (4.20) 


Simulation 


Mean queue length 


1.475 


6.359 


6.413 


Abandonment fraction 


0.05863 


0.05517 


0.05512 


P[jV(oo) > 480] 


0.8663 


0.8929 


0.8881 


P[JV(co) > 500] 


0.3192 


0.4822 


0.4720 


P[AT(oo) > 520] 


9.274 x 10" 5 


0.1074 


0.1050 


P[iV(oo) > 550] 


-4.488 x 10" 9 


0.006616 


0.006248 



(b) p = 1.045 and n = 500 
Table 5 

Performance measures of the M/H2/11 + H2 queue. 



times follow a two-phase hyperexponential distribution that has 

p= (0.9,0.1)' and v = (1,200)'. 

In other words, there are two types of patience times. Ninety percent of 
patience times are exponentially distributed with mean one and ten percent 
are exponentially distributed with mean 0.005. We take the same hyperex- 
ponential service time distribution as in Sections 6.2 and 6.3. 

The hazard rate function of the hyperexponential patience time distribu- 
tion is 

hit) = P^eM-ht) +p 2 heM-V2t) for f > Q _ 
pi exp(—vit) + p2 exp(-i> 2 t) 

The drift coefficient of X in (4.20) is also given by (6.4) where 



log (pi exp ( —i>iz) + P2 exp ( 7-^2% ) ) for z > 0. 



in 

In this example, we have 

h(0) =pih+P2i>2 = 20.9 and /i (1) (0) = -pip 2 (h - V2? = -3564.1. 



Thus, £q = and the hazard rate function has a steep slope near the ori- 
gin. Since the zeroth degree Taylor's approximation in (5.9) may bring on 
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too much error, the reference density exploiting the lowest-order nonzero 
derivative at the origin could be erroneous. 

To choose an appropriate reference density, an auxiliary queue is used 
again. As in Section 5.3, the auxiliary queue is an MjH%ln + M queue 
that shares the same arrivals and service times with the M/H^/n + H2 
queue. Let a > be the rate of the exponential patience time distribution. 
We take a = v\ A v<i so that the patience times in the auxiliary queue 
are all of the type with the longer mean. If the queue lengths are equal, 
the abandonment rate in the auxiliary queue must be lower than that in 
the original queue. Therefore, the queue length decays slower in the former 
queue and the reference density for model (4.14) of the auxiliary queue 
should work. This observation leads to a reference density that follows (5.6) 
and (5.14), but in this example, we take a = i>\ A z>2 and solve (5.11) to find 
<7o- 

Two scenarios with n = 50 and 500 servers are investigated. The respec- 
tive arrival rates are A = 57.071 and 522.36. Thus, p = 1.141 and 1.045 and 
both scenarios have f3 = —1. By solving (5.11), we have qo = 0.165 for the 
first scenario and go = 0.0059 for the second scenario. The reference density 
follows (5.6) and (5.14) with a = v\ = 1. With eq = 10 -7 , the truncation 
rectangle is K = [—7,9] x [—7,9], partitioned into 0.5 x 0.5 finite elements. 
The performance estimates obtained by the diffusion model (4.20) are com- 
pared with the simulation results in Table 5. The performance estimates are 
still quite accurate. 

We also put the performance estimates produced by the diffusion model 
(4.14) in this table. For this model, the reference density follows (5.6) and 
(5.8) with a = h(0) = 20.9. In the proposed algorithm, the mesh for model 
(4.20) is used again. Because in this example, using the patience time density 
at zero solely cannot capture the behavior of the abandonment process, 
model (4.14) fails to produce proper performance estimates. 

7. Implementation issues. The proposed algorithm was implemented 
using the C++ programming language. The package was tested on both 
Linux and Windows platforms. In this section, we discuss several important 
issues arising from the implementation. They are crucial for using the algo- 
rithm to solve practical problems. To demonstrate these issues, the second 
scenario with n = 500 servers in Section 6.1 is investigated throughout this 
section. The diffusion model (4.14) is used to approximate the MjHiju+M 
queue. 

7.1. Influence of the reference density. The reference density plays a key 
role in the algorithm. If the function r does not satisfy (3.1), the sequence 
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400 500 600 700 800 400 500 600 700 800 

the number of customers in system the number of customers in system 



(a) K = [-7, 10] x [-7, 10] (b) K = [-7, 32] x [-7, 32] 

Fig 5: The output of the proposed algorithm with the "naive" reference 
density. 

of spaces {H^ : k G N} may not converge to H in L 2 (M rf ,r) and the output 
of the algorithm may significantly deviate from the exact stationary density. 
To demonstrate this issue, let us consider a "naive" reference density. 

To produce a "naive" reference density, we consider a queue that has the 
same arrival process and patience time distribution as the M/H2/n + M 
queue. This new queue has an exponential service time distribution and its 
mean service time is equal to that of the M/H2/TI + M queue. For this 
M/M/n + M queue, the diffusion model (4.14) is a one-dimensional piece- 
wise OU process whose stationary density is given by (5.5). The "naive" 
reference density is a product reference density in (5.6) with each rj being 
the stationary density in (5.5). In other words, the "naive" reference density 
is obtained by pretending the service time distribution to be exponential. 

Let us apply the "naive" reference density. With eo = 1CP 7 , the trunca- 
tion rectangle is set to be K = [—7,10] x [—7,10] and is partitioned into 
0.5 x 0.5 finite elements. As shown in Figure 5a, the output of the pro- 
posed algorithm noticeably deviates from the exact stationary distribution. 
To further confirm that the "naive" reference density cannot work, we next 
test the truncation rectangle K = [—7,32] x [—7,32], which is used in Sec- 
tion 6.1 along with the proposed reference density. In this case, the matrix 
A in (3.12) is close to singular and its condition number is 3.52 x 10 190 . 
Figure 5b manifests the severe error in the algorithm output. 

Recall that in this example, the hyperexponential service time distribution 
has c 2 . = 24. Comparing (5.5) with (5.8), we can tell that the decay rate 
of the "naive" reference density is much larger than that of the proposed 
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-proposed algorithm 
- matrix-analytic 



proposed algorithm 

matrix-analytic 



the number of customers in system 



the number of customers in system 



(a) 1.0 x 1.0 finite elements (b) 0.25 x 0.25 finite elements 

Fig 6: The output of the proposed algorithm with different meshes. 





0.5 x 0.5 


0.25 x 0.25 


Matrix-analytic 


Mean queue length 


54.17 


54.17 


54.05 


Abandonment fraction 


0.05181 


0.05182 


0.05173 


P[JV(oo) > 470] 


0.9701 


0.9702 


0.9694 


P[N(oo) > 500] 


0.6838 


0.6835 


0.6818 


P[N(oo) > 600] 


0.2244 


0.2241 


0.2229 


P[AT(oo) > 750] 


0.008233 


0.008246 


0.006395 



Table 6 



The output of the proposed algorithm using different meshes. 



reference density. If Conjecture 3 is true, one can expect that the "naive" 
reference density decays much faster than the stationary density and the 
second condition in (3.1) may not hold. In this case, the ratio function q is 
no longer in L 2 (M. d ,r) and consequently, the algorithm fails to produce any 
adequate estimate of the ratio function. 

7.2. Mesh selection. When both the reference density and the truncation 
hypercube are fixed, using a finer mesh may produce smaller approximation 
error. However, a finer mesh yields more basis functions, which in turn lead 
to a larger condition number for the matrix A in (3.12). If the condition 
number of A is too large, the round-off error in solving (3.12) becomes 
considerable. So a finer mesh does not necessarily yield a more accurate 
output. 

Let us test different meshes for the second scenario in Section 6.1. We keep 
the same settings for the algorithm except the size of finite elements. The 
output with 1.0 x 1.0 finite elements is plotted in Figure 6a. With this mesh, 
the algorithm does not perform well at the intervals where the stationary 
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400 500 600 700 800 400 500 600 700 800 

the number of customers in system the number of customers in system 



(a) 0.25 x 0.25 finite elements (b) 0.125 x 0.125 finite elements 

Fig 7: The output of the proposed algorithm with the "naive" reference 
density and different meshes. 

density varies quickly. We need a finer mesh to improve the accuracy. In this 
case, the condition number of A is 5.70 x 10 20 . Recall that to produce the 
curve in Figure lb, we use a mesh consisting of 0.5 x 0.5 finite elements. With 
this mesh, the condition number of A is 1.15 x 10 23 . When the element size 
is further reduced to 0.25 x 0.25, the condition number of A grows to 7.13 x 
10 27 . As illustrated in Figure 6b, the output of the algorithm fits the exact 
stationary distribution well. When we compare Figures lb and 6b, however, 
there is barely any difference noticeable between the algorithm outputs. To 
confirm that this mesh is not superior to the one with 0.5x0.5 finite elements, 
we list several performance estimates in Table 6. In this table, the results 
in Table lb are duplicated for comparison purposes. The difference between 
the algorithm outputs using these two meshes is negligible. Considering the 
modeling error of the diffusion model, we can assert that using 0.5 x 0.5 
finite elements is sufficient to produce an accurate approximation for this 
queue. 

Given an appropriate reference density and the associated truncation hy- 
percube, the above discussion has demonstrated an approach to selecting a 
mesh. Beginning with two meshes, with one finer than the other, we com- 
pare the algorithm outputs using these two meshes. If obvious difference is 
observed, the coarser mesh should be discarded and a further finer mesh is 
explored. Continue this procedure until the difference between the outputs 
of two meshes are negligible. Then, the coarser one of the remaining two is 
selected as an appropriate mesh. 

We would also demonstrate that with an improper reference density, a 
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m = 4 


m = 8 


m — 16 


Matrix-analytic 


Mean queue length 


54.17 


54.17 


54.17 


54.05 


Abandonment fraction 


0.05181 


0.05181 


0.05181 


0.05173 


P[Ar(oo) > 470] 


0.9701 


0.9701 


0.9701 


0.9694 


P[JV(oo) > 500] 


0.6833 


0.6838 


0.6839 


0.6818 


P[JV(oo) > 600] 


0.2245 


0.2244 


0.2244 


0.2229 


P[AT(oo) > 750] 


0.008235 


0.008233 


0.008232 


0.006395 



Table 7 

The output of the proposed algorithm with different quadrature orders. 



finer mesh cannot make the algorithm yield an adequate output. Let us go 
back to the example in Section 7.1 with the "naive" reference density. We 
set the truncation rectangle to be K = [—7,32] x [—7,32] and the size of 
finite elements to be 0.25 x 0.25. The output is shown in Figure 7a. Although 
the curve appears smoother than the one in Figure 5b with 0.5 x 0.5 finite 
elements, the output still fails to capture the exact stationary distribution. 
This time, the condition number of A is 3.91 x 10 195 . There is no doubt that 
such an ill-conditioned matrix will bring about a huge round-off error in 
solving (3.12). A mesh with 0.125 x 0.125 finite elements is also investigated 
and the algorithm output is plotted in Figure 7b. The condition number of 
A increases to 6.35 x 10 198 and the algorithm misses the target as well. 

7.3. Gauss-Legendre quadrature. Before solving the linear system (3.12), 
we must generate the matrix A and the vector v whose entries are given by 
(3.13). We follow a Gauss-Legendre quadrature rule to compute the integral 
for each entry. The integral is taken over a two-dimensional rectangle and 
the quadrature rule evaluates the integrand at m points in each dimension. 
The results are more accurate when a larger m is used. In Section 6, we take 
m = 8 in the numerical examples. Here, we briefly discuss the impact of the 
order m. 

Several performance estimates are listed in Table 7. We keep the same 
settings for the algorithm except the quadrature order in each dimension. 
For the convenience of comparison, the results in Table lb are duplicated 
in this table. Clearly, the Gauss-Legendre quadrature of order m > 4 is 
sufficiently accurate for our purposes. 

7.4. Computational complexity. Let d, the dimension of the diffusion 
model, be fixed. The size of A is mc X mc where mc is the dimension 
of the functional space C given by (3.18). The matrix A is sparse. There are 
at most 6 d nonzero entries in each row or column. Hence, it takes 0(mc) 
arithmetic operations to construct A. We may apply Gaussian elimination 
to solve the linear system (3.12). When the basis functions are properly or- 
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1.0 x 1.0 


0.5 x 0.5 


0.25 x 0.25 


0.125 x 0.125 


Dimension mc 


5776 


23716 


96100 


386884 


Constructing A and v 


6.63 


27.3 


109 


455 


Solving (3.12) 


0.0780 


0.359 


2.29 


18.2 



Table 8 



Computation time (in seconds) of the proposed algorithm using different meshes. 

dered, the nonzero entries of A are confined to a diagonally bordered band of 
width 0(m^ 1 ^ d ). Hence, solving (3.12) requires 0(m^, d *" d ) operations 
as mc — > oo. 

The computation time (measured by seconds) for various meshes can be 
found in Table 8, where we list both the time for constructing A and v 
and the time for solving (3.12). When computing A and v, we follow a 
Gauss-Legendre quadrature rule with m = 8 points in each dimension. The 
truncation rectangle is set to be K = [—7,32] x [—7,32]. Each mesh is 
obtained by setting the size of finite elements. The dimension mc increases 
by around four times as the width of each finite element is reduced by half. 
The proposed algorithm is tested on a laptop with a 2.66GHz Intel Core 2 
Duo processor and eight gigabytes memory. Both A and v are produced by 
our C++ package. The linear system (3.12) is solved by Matlab. These two 
parts are connected via a MEX interface that comes with Matlab. 

8. Concluding remarks. In this paper, we proposed two approximate 
models for many-server queues with customer abandonment. Both these 
models are diffusion processes and they differ in how the abandonment pro- 
cess is approximated. A finite element algorithm was proposed for computing 
the stationary distribution of each model. The essential part of the algorithm 
is a reference density that controls the convergence of the algorithm. To con- 
struct a reference density, we conjectured that the limit queue length process 
has a certain Gaussian tail. Using this conjecture, we proposed a systematic 
approach to choosing a reference density. With the proposed reference den- 
sity, the output of the algorithm is stable and accurate. Numerical examples 
indicate that the diffusion models are good approximations for many-server 
queues. 

Assume that the stationary density g is twice differentiable in M. d and 
vanishes at infinity. Using the basic adjoint relationship (2.7) and applying 
integration by parts twice, we have 

g*g(x) = for all x£R d 

where Q* is the adjoint operator of the generator Q. Fix a finite domain 
Kcl' ! large enough. One can solve the stationary density g by the Dirichlet 
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problem 

{G*g(x) = for x in the interior of K, 
g(x) = for x on the boundary of K. 

Such a Dirichlet problem can be solved via a finite difference algorithm. 
Alternatively, for each test function /, one may apply integration by parts 
once to the basic adjoint relationship to obtain an equation that involves 
the first derivatives of g and the first derivatives of /. From this weak for- 
mulation, fixing a large enough finite domain K and assuming that g is 
zero on the boundary of K, one may apply a standard Galerkin finite ele- 
ment method to compute the stationary density g on K. See, e.g., Kovalov, 
Linetsky and Marcozzi (2007). Both the finite difference algorithm and the 
Galerkin method do not use a reference density. A future research topic is 
to compare the efficiency and accuracy of these two algorithms with the 
proposed algorithm in this paper. 

The dimension of the functional space C in Section 3.3 grows exponen- 
tially in d, the dimension of the diffusion model. As a consequence, both 
the computation time and the memory usage increases exponentially in d. 
When d is not small, the curse of dimensionality is a serious challenge for the 
proposed algorithm as well as any other algorithms. To reduce the dimen- 
sion of C, one possible approach is to investigate a reference density that 
potentially shares more common features with the stationary density. Such 
a reference density may enable us to compute the stationary density with a 
moderate number of basis functions when d is not small. Another possible 
direction to reduce the computational complexity of the algorithm is to in- 
vestigate a low-rank matrix approximation for the linear system (3.12). The 
technique of random sampling may be explored. See Kannan and Vempala 
(2010) for more details. 

APPENDIX 

Proof of Proposition 2. Recall that K is the compact support of C 
and the basis functions of C are given by (3.17). We use Cq(K) to denote the 
set of real-valued functions on a neighborhood of K that are continuously 
differentiable and have compact support in K. Clearly, C C Cq(K). For any 
/, / G Cq(K), we define an inner product by 



3=1 

rl,2 



and let W ' (K) be the closure of Q}(iT) in the norm induced by this inner 
product. Then, W ' (K) is a Hilbert space and C C Wq' z (K). 
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Proof of Proposition 2. Since Q is a linear operator, it suffices to 
show that for any fo G C, we must have fo = if £/_fo = in L 2 (R d ,r). 
The uniform elliptic operator Q can be written into the divergence form as 
in (8.1) of Gilbarg and Trudinger (2001), i.e., 

e/W = fi 3 (# + -ff d{^{ x )dj( x) /d X] ) 

3=1 3 3=1 1=1 

for each / G (IR '), where 

^) = W-2g-^- 

Let U C l d be a connected open set that is bounded and contains K. Since 
r > and Q fo is continuous in the interior of each finite element, we must 
have Gfo = in K except on the boundaries of certain finite elements where 
Qfo is not defined. Hence, Qfo = in U in the weak sense (see (8.2) of Gilbarg 
and Trudinger (2001)). Note that b, E, and the partial derivatives of S are 
all continuous, so both b and S are bounded in U. Because fo G W Q ' (K), 
it follows from Corollary 8.2 of Gilbarg and Trudinger (2001) that fo = in 
K, and thus / = in R d . □ 

Proof of Proposition 3. Given a compact set K C let C^(K) 
be the set of real-valued functions on a neighborhood of K that are twice 
continuously differentiable with bounded first and second derivatives in K. 
For each / G C^(K), define a norm ||-||#2(#) by 



Because both b and S are bounded in K, there exists kq{K) > such that 
(A.l) j^gf(x)) 2 r(x)dx < K (K) \\f\\ 2 mK) for all / G C^K). 



Let C 2 (K) be the closure of C 2 (K) in the above norm. A standard procedure 
can be used to define the first-order and the second-order derivatives for each 
/ G C 2 {K). Then, the operator Q can be extended to C 2 {K) and inequality 



(A.l) holds for all / G C?(K). 



Proof of Proposition 3. It suffices to prove that for any / G C 2 (R d ), 
there exists a sequence of functions {ip^ G : k G N} such that 

\\QVk - Qfo\\ ->■ as k oo. 
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Fix e > 0. Because f Mr as k — > oo, by (3.2) and the Cauchy-Schwartz 
inequality, there exists a G N such that 



(A.2) / (<?/ (x)) 2 r(x)d:r 



e 2 
< 2 



Consider the finite hypercube K a . By (A.l), there exists Ko(K a ) > such 
that 

(A.3) f (Gf{x)) 2 r{x) dx < K (K a ) ||/|| W a) for all / G (7 2 (K a ). 

A polynomial can be used to approximate /o on K a . By Proposition 7.1 in 
the appendix of Ethier and Kurtz (1986), there exists a polynomial f p such 
that 

For the lattice mesh Ajt, let A ai fc be the set of its nodes in the interior of 
K a . For any k > a, let be a function in such that <fk(z) = for all 
x€R d \K a and 

= / P (x) and *g^> = d JM 

for j = 1, . . . , d and all x G A Qj fc. Clearly, G C 2 (K a ). Because the sequence 
of meshes {A^ : k G N} is regularly refined, there exists a constant K\ > 
such that 77a*. < ^1 for all k > a. Using the interpolation error estimate in 
Theorem 6.6 of Oden and Reddy (1976), we have 

2 ( f \ 1/2 2 

llWfc - fp\\H2(K ) ^ K i K 2K3( / r(x)dx) |A fc | , 

where K2 > is a constant independent of A& and / p , and 
d 4 f P (x) 



K3 = sup 



dx™ 1 • • • dx™ d 



x G K a ; mi + • • • + m<2 = 4 > < 00. 



Hence, there exists <5o > such that 

whenever |A^| < 5$. In this case, 

\\<Pk ~ fo\\ H \K a ) < II V* - fv\\ H *{K a ) + ll/p " /ollff2(tf a ) < ~ /|=^ ' 
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By (A.3) 



(A.4) / (g<p k (x) - Gfo(x)) 2 r(x) dx < K (K a ) \\ip k - f 



JK a 



H H 2 {K a ) 




It follows from (A. 2) and (A.4) that \\Gfk — GfoW < e whenever k > a and 



The authors would like to thank Ton Dieker and Vadim Linetsky for their 
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